######################################################## 
## /projects/bsi/hematology/s106284.gwas/rpgm/beagle_filenum.R
## Last Modified: Friday, 10 September 2010 11:30 AM CDT
######################################################## 
##  Program Name  :/projects/bsi/hematology/s106284.gwas/rpgm/beagle_filenum.R
##  Study Title   : GENEVA-Heit
##  Programmer    : MEM
##  Investigator  : de Andrade
##  Statistician  : Armasu,Petterson
##  Creation eDate: Monday, 08 March 2010 11:26 AM CST
##  Study Number  : s106284
##  Function      : For BEAGLE Chr X imputation, make phased files
## Break up into males/females, clean up files from HapMap
######################################################## 
hapdirX <- "/data1/bsi/hematology/s106284.gwas/impute/hapmap_chrX"

setwd(hapdirX)
reffiles <- list.files(pattern="phased$")
reffiles

for (file in reffiles) {

  filename <- gsub(".*chr23_", "", gsub(".phased", "", file))

  temp <- read.table(file, header=T, as.is=T)

  assign(filename, temp)

  rm(temp)
  rm(filename)
}

##Get HapMap populations
pop <- read.table("/data1/bsi/hypertension/s104892.genoa/plinkfiles/hapmap/hapmap_subj.txt", header=T, as.is=T)

ls()
ceu.unr[1:5,1:5]
ceu.trios[1:5,1:5]
ceu.duos[1:5,1:5]
tsi.unr[1:5,1:5]

dim(ceu.unr)
dim(ceu.trios)
dim(ceu.duos)
dim(tsi.unr)

all <- merge(ceu.unr, ceu.trios, by=c("rsID", "position_b36"),all=T)
dim(all)
all <- merge(all, ceu.duos, by=c("rsID", "position_b36"),all=T)
dim(all)
all <- merge(all, tsi.unr, by=c("rsID", "position_b36"),all=T)
dim(all)

##Get ids in phased data
refids <- names(all)[-(1:2)]
refids[1:5]
##These are unobserved parent haps per readme file, drop per Sooraj. 
refids[-grep("^NA[0-9]{5}_(A|B)$", refids)]

##Keep ids
reftab <- table(gsub("_(A|B)", "", refids[grep("NA[0-9]{5}_(A|B)", refids)]))
table(reftab)
keepids <- paste(rep(names(reftab)[reftab==2],each=2), rep(c("_A", "_B"), length(reftab[reftab==2])),sep="")

all <- all[,c("rsID", "position_b36", keepids)]
dim(all)

##Are all members of reference sample founders?
##NA10850 has momid ne 0, but mother is not in dataset
table(pop$dad[pop$IID %in% names(reftab)[reftab==2]])
table(pop$mom[pop$IID %in% names(reftab)[reftab==2]])
pop[pop$mom=="NA12058",]
pop[pop$IID=="NA12058",]
pop[pop$FID==1344,]
reftab[grep("(NA12058|NA10850|NA10851|NA12056|NA12057)", names(reftab))]

##Put in BEAGLE format
all <- all[order(all$"position_b36"),]
all <- cbind("I"=rep("M", nrow(all)), all, stringsAsFactors=F)
names(all)[names(all)=="rsID"] <- "id"

keepids.f <- keepids[gsub("_(A|B)", "" ,keepids) %in% pop$IID[pop$sex==2]]
keepids.m <- keepids[gsub("_(A|B)", "" ,keepids) %in% pop$IID[pop$sex==1]]
length(keepids.f)
length(keepids.m)

#Separate out pseudo-autosomal region
#Region 1: ChrX:1-2709520 (ChrY:1-2709520)
#Region 2: ChrX:154584237-154913754 (ChrY: 57443437-57772954)
pseudmark <- all[(all$"position_b36"<=2709520)|((all$"position_b36">=154584237)&(all$"position_b36"<=154913754)),"id"]
length(pseudmark)

all.f <- all[,c("I", "id", "position_b36", keepids.f)]
all.m <- all[,c("I", "id", "position_b36", keepids.m)]

##Get rid of second hap for males
allmiss <- apply(all.m!="-",2,sum)==0
length(allmiss)
dim(all.m)
all.m2 <- all.m
all.m <- all.m[,(1:ncol(all.m))[!allmiss]]
all.m[1:10,1:10]
keepids.m <- names(all.m)[grep("NA[0-9]{5}_(A|B)", names(all.m))]
sum(gsub("_(A|B)", "", names(all.m2)[which(allmiss)])==gsub("_(A|B)", "", keepids.m))
length(keepids.m)
all.m2[,which(allmiss)] <- all.m2[,keepids.m]

##Get unique alleles for each marker
all.f$allA <- apply(all.f[,keepids.f]=="A",1,sum)>0
all.f$allC <- apply(all.f[,keepids.f]=="C",1,sum)>0
all.f$allG <- apply(all.f[,keepids.f]=="G",1,sum)>0
all.f$allT <- apply(all.f[,keepids.f]=="T",1,sum)>0
all.f$markers <- gsub("[[:space:]][[:space:]]*", " ",paste(c("", "A")[all.f$allA+1],c("", "C")[all.f$allC+1],c("", "G")[all.f$allG+1],c("", "T")[all.f$allT+1]))
all.f$markers <- gsub("[[:space:]]$", "", all.f$markers)
all.f$markers <- gsub("^[[:space:]]", "", all.f$markers)
table(all.f$markers)

all.m$allA <- apply(all.m[,keepids.m]=="A",1,sum)>0
all.m$allC <- apply(all.m[,keepids.m]=="C",1,sum)>0
all.m$allG <- apply(all.m[,keepids.m]=="G",1,sum)>0
all.m$allT <- apply(all.m[,keepids.m]=="T",1,sum)>0
all.m$markers <- gsub("[[:space:]][[:space:]]*", " ",paste(c("", "A")[all.m$allA+1],c("", "C")[all.m$allC+1],c("", "G")[all.m$allG+1],c("", "T")[all.m$allT+1]))
all.m$markers <- gsub("[[:space:]]$", "", all.m$markers)
all.m$markers <- gsub("^[[:space:]]", "", all.m$markers)
table(all.m$markers)

table(all.m$markers, all.f$markers)

markdat <- all.f[,c("id", "position_b36", "markers")]
markdat$markers.m <- all.m$markers[match(markdat$id, all.m$id)]
table(nchar(markdat$markers))
markdat$markers[nchar(markdat$markers)==1] <- markdat$markers.m[nchar(markdat$markers)==1]
table(nchar(markdat$markers))
markdat <- markdat[nchar(markdat$markers)>1,]
dim(markdat)

##Write out marker data
setwd("/data1/bsi/bioinf_int/s106381.borawork/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/beagle_mainrun/workflow/files/reference/hapmapfiles/xchr/test")
write.table(markdat[,c("id", "position_b36", "markers")], row.names=F, col.names=F, quote=F, "markers.txt")
system("gzip markers.txt")

##Subset to non-monomorphic markers
all.f <- all.f[all.f$id %in% markdat$id,c("I", "id", keepids.f)]
all.m <- all.m[all.m$id %in% markdat$id,c("I", "id", keepids.m)]
all.m2 <- all.m2[all.m2$id %in% markdat$id,c(1,2,grep("NA[0-9]{5}_(A|B)", names(all.m2)))]

dim(all.f)
dim(all.m)
dim(all.m2)

all.f <- rbind(gsub("_(A|B)$", "", names(all.f)), all.f)
all.m <- rbind(gsub("_(A|B)$", "", names(all.m)), all.m)
all.m2 <- rbind(gsub("_(A|B)$", "", names(all.m2)), all.m2)

all.f[1:5,1:5]
all.m[1:5,1:5]

write.table(all.f,"phasedfemale.bgl", row.names=F, col.names=F, quote=F)
system("gzip phasedfemale.bgl")

write.table(all.m,"phasedmale.bgl", row.names=F, col.names=F, quote=F)
system("gzip phasedmale.bgl")

write.table(all.m2,"phasedmaledip.bgl", row.names=F, col.names=F, quote=F)
system("gzip phasedmaledip.bgl")

##Read in gwas subjects, separated into males and females
keepsubj <- read.table("/data1/bsi/hematology/s106284.gwas/txt/GENEVA_Eigensoft_6SD_keepsubj.txt", header=F, as.is=T)
fam <- read.table("/data1/bsi/hematology/s106284.gwas/Plink/heitassocgwa.fam", header=F, as.is=T)

keepsubj.m <- keepsubj[keepsubj$V1 %in% fam$V1[fam$V5==1],]
keepsubj.f <- keepsubj[keepsubj$V1 %in% fam$V1[fam$V5==2],]
dim(keepsubj.m)
dim(keepsubj.f)

write.table(keepsubj.m, row.names=F, quote=F, col.names=F, "gwas_keepsubj_males.txt")
write.table(keepsubj.f, row.names=F, quote=F, col.names=F, "gwas_keepsubj_females.txt")

##Check strand
#map <- read.table("/data1/bsi/hematology/s106284.gwas/Plink/heitassocgwa.bim", header=F, as.is=T)
map <- read.table("/data1/bsi/hematology/s106284.gwas/Plink/merged/mergedassocgwa.bim", header=F, as.is=T)
names(map) <- c("chr", "rsid", "cm", "bp", "a1", "a2")
map <- map[map$rsid %in% markdat$id,]
dim(map)
map[1:5,]
map$alleles <- apply(map[,c("a1", "a2")], 1, function(x) {paste(sort(x)[1], sort(x)[2], sep=" ")})

map$bgl <- markdat$markers[match(map$rsid, markdat$id)]
table(map$alleles, map$bgl)

check <- map[(map$alleles!=map$bgl)|(map$alleles %in% c("A T", "C G")),]
check <- check[check$alleles!="0 0",]
check

flip <- check[!(check$alleles %in% c("A T", "C G")),]
dim(flip)
table(flip$alleles, flip$bgl)
err <- flip[!((flip$alleles=="A C"&flip$bgl=="G T")|(flip$alleles=="A G"&flip$bgl=="C T")|(flip$alleles=="C T"&flip$bgl=="A G")|(flip$alleles=="G T"&flip$bgl=="A C")),]
flip <- flip[((flip$alleles=="A C"&flip$bgl=="G T")|(flip$alleles=="A G"&flip$bgl=="C T")|(flip$alleles=="C T"&flip$bgl=="A G")|(flip$alleles=="G T"&flip$bgl=="A C")),]
check <- check[(check$alleles %in% c("A T", "C G")),]
check <- rbind(check, err)
dim(flip)
dim(err)
dim(check)

##Get frq of CG, AT, 0 SNPs
#freq <- read.table("/data1/bsi/hematology/s106284.gwas/output/all.frq", header=T, as.is=T)
freq <- read.table("/data1/bsi/hematology/s106284.gwas/output/merged.frq", header=T, as.is=T)
freq <- freq[freq$SNP %in% check$rsid,]
dim(freq)
freq

checksnp <- function(snp) {
  hapfrq <- table(as.character(all[all$id==snp,keepids]))
  hapfrq <- hapfrq[names(hapfrq)!="-"]
  hapfrq <- hapfrq/sum(hapfrq)

  plinkfrq <- freq[freq$SNP==snp,c("SNP", "A1", "A2", "MAF")]

  return(list("hap"=hapfrq, "data"=plinkfrq))
}

check
checksnp("rs5948755") ##okay
checksnp("rs4830628") ##okay
checksnp("rs6633383") ##okay
checksnp("rs12009538")  ##flip
checksnp("rs6568006") ##Drop? Looks weird, A should be minor allele but is monomorphic A
checksnp("rs1858934") ##Drop, AT SNP, MAF close to .5, can't tell if right strand or not
checksnp("rs6626075") ##okay
checksnp("rs4898352") ##flip

setwd(hapdirX)
write.table(c("rs6568006", "rs1858934"), "drop_impute.txt", row.names=F, col.names=F, quote=F)
#write.table(c(flip$rsid, "rs12009538"), "flip_phase3.txt", row.names=F, col.names=F, quote=F)
write.table(c(flip$rsid, "rs12009538", "rs4898352"), "flip_phase3_merged.txt", row.names=F, col.names=F, quote=F)
